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It is usually assumed that relic neutrinos possess a Fermi-Dirac distribution, acquired during ther- 
mal equilibrium in the Early Universe. However, various mechanisms could introduce strong dis- 
tortions in this distribution. We perform a Bayesian likelihood analysis including the first moments 
of the three active neutrino distributions as free parameters, and show that current cosmological 
observations of light element abundances, Cosmic Microwave Background (CMB) anisotropics and 
Large Scale Structures (LSS) are compatible with very large deviations from the standard picture. 
We also calculate the bounds on non-thermal distortions which can be expected from future observa- 
tions, and stress that CMB and LSS data alone will not be sensitive enough in order to distinguish 
between non-thermal distortions in the neutrino sector and extra relativistic degrees of freedom. 
This degeneracy could be removed by additional constraints from primordial nucleosynthesis or 
independent neutrino mass scale measurements. 
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I. INTRODUCTION 

The large amount of observations of the Cosmic Mi- 
crowave Background (CMB) which have been accumu- 
lated through forty years, since its discovery till the most 
recent detailed CMB map provided by the Wilkinson Mi- 
crowave Anisotropy Probe (WMAP) Collaboration 
shows that relic photons are distributed according to a 
black body distribution to an incredibly high accuracy. 
This fact is indeed one of the main cornerstones of the 
hot Big Bang model, strongly supporting the idea that 
our Universe expanded starting from high density and 
temperature conditions. 

It is well known that in this framework we expect the 
Universe to be also filled by a large amount of neutri- 
nos, with present densities of the order of 10 2 cm~ 3 per 
flavor. These relic neutrinos decoupled from the electro- 
magnetic plasma quite early in time, when weak interac- 
tion rates became slower than the Hubble rate for (pho- 
ton) temperatures in the range 2-^4 McV, just before 
Big Bang Nucleosynthesis (BBN) took place. Unfortu- 
nately the fact that this Cosmic Neutrino Background 
(CzaB) has today a very small kinetic energy, of the order 
of 10~ 4 eV, and that neutrinos only interact via weak 
interactions, prevents us from any possible direct detec- 
tion of this background on the Earth (see e.g. 0, 0] for 



recent reviews). Nevertheless, there are several indirect 
ways to constrain the CiaB by looking at cosmological 
observables which are influenced either by the fact that 
neutrinos contribute to the Universe expansion rate at 
all stages, or also via their interactions with the electro- 
magnetic plasma and baryons before their decoupling. 
In this respect the most sensitive probe is represented 
by the values of the light nuclide abundances produced 
during BBN. Actually, the final yields of Deuterium, 7 Li 
and in particular of He strongly depend on the num- 
ber of neutrino species as well as on their distribution in 
phase space at about 1 MeV when the neutron to pro- 
ton density ratio freezes. In fact, since neutrinos were 
in chemical equilibrium with the electromagnetic plasma 
till this epoch we know by equilibrium thermodynamics 
that they were distributed according to a Fermi-Dirac 
function, yet BBN can constrain exotic features like the 
value of their chemical potential 0,|EIE0- Furthermore, 
a detailed analysis carried out by several authors 0, 
shows that neutrinos benefit to a small extent of the en- 
tropy released by electron-positron pairs at their annihi- 
lation stage, and that this phenomenon reheat neutrinos 
in a non-thermal way, introducing a small distortion of 
the order of few percent which grows with the neutrino 
comoving momentum. 

Summarizing, neutrinos arc likely to emerge from the 
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BBN epoch as decoupled species with an equilibrium dis- 
tribution, a smaller temperature with respect to photons 
by the approximate factor (4/11) 1 / 3 , possibly a small 
chemical potential to temperature ratio |£| < 0.1 1 and 
finally, tiny momentum dependent non-thermal features. 

Once decoupled, neutrinos affect all key cosmological 
observables which are governed by later stages of the evo- 
lution of the Universe only via their coupling with gravity. 
A well known example is their contribution to the total 
relativistic energy density, which affects the value of the 
matter-radiation equality point, which in turn influences 
the CMB anisotropy spectrum, in particular around the 
scale of the first acoustic peak. Similarly, their number 
density and their masses are key parameters in the small 
scale suppression of the power spectrum of Large Scale 
Structures (LSS) , smaller than l nr , the horizon when neu- 
trinos become non-relativistic [lOj 
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where tu m = fl m h 2 , with Cl m the fraction of the crit- 
ical density due to matter and h the Hubble constant 
in units of 100 km s" 1 Mpc" 1 . Several authors have 
investigated in details these issues assuming a stan- 
dard Fermi-Dirac distribution for neutrinos, see e.g. 
[1 El El El El El Indications for primordial 
anisotropies in the neutrino distribution function have 
been pointed out in |l7j . In this paper we consider the 
possibility that the relic neutrino distribution in phase 
space may be sensibly different and address the follow- 
ing point: how present (and future) cosmological obser- 
vations can prove that indeed neutrinos are thermally 
distributed? The implications of a relic neutrino spec- 
trum significantly different from a Fermi-Dirac distribu- 
tion have been also considered by several authors. For 
instance, neutrinos could violate the Pauli exclusion prin- 
ciple and obey Bose-Einstein statistics with important 
cosmological and astrophysical implications El El kSI • 
Non-equilibrium neutrino spectra are also produced in 
low-reheating scenarios I22I |23| , from the decays of 
massive neutrinos into relativistic products (see for in- 
stance pil [25|]) or as the result of active-sterile neu- 
trino oscillations after decoupling (see [2|J and references 
therein) . A full list of references can be found in j23| ■ 

In order to answer the question raised above, we con- 
sider in detail (Section [nj a simple model where the 
out of equilibrium decay of a light scalar produces non- 
thermal features in the neutrino spectra, although our 
results are quite general and can be also applied to dif- 
ferent scenarios. We describe the current bounds on this 



1 We point out that the bound on £ is shared by all neutrino species 
due to flavor oscillations and that a much looser bound is 
obtained if we allow for extra relativistic degrees of freedom con- 
tributing to the Hubble rate during the BBN epoch (see e.g. 
0). 



model in Section HTH both from BBN and CMB+LSS. 
The expected sensitivity of future cosmological data can 
be forecast if a fiducial cosmological model is assumed, as 
described in Section llVl Finally we report our concluding 
remarks in Section Ivl 



II. THE PHASE-SPACE DISTRIBUTION OF 
RELIC NEUTRINOS 

As we already said, till temperatures of the order of 
MeV neutrinos are in chemical equilibrium with the elec- 
tromagnetic plasma, so that they keep an equilibrium 
distribution. After weak interaction decoupling and in 
absence of any other interaction process, the collision- 
less kinetic equations in a Lemaitre- Friedman- Robertson- 
Walker universe state that, assuming for simplicity in- 
stantaneous decoupling, the neutrino and antincutrino 
distribution in phase-space is a Fermi-Dirac function. 
Therefore, in terms of the comoving momentum y = ka, 
with a the scale factor and neutrino temperature T v , one 
has 
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for each mass eigenstate. We assume zero chemical po- 
tential in the following. Regardless of the specific case at 
hand, we can specify the neutrino distribution df a (y) by 
the set of moments Qa 
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where T is the photon temperature, and we have de- 
fined the in terms of the standard value of the neu- 
trino temperature in the instantaneous decoupling limit 
T„ = (4/ll) 1 / 3 T. We assume in the following that neu- 
trino distribution decays at large comoving momentum 
as exp(— y), so that the distribution admits moments of 
all orders. Actually, this is not a particularly severe con- 
straint, since we expect that at very high y the shape of 
the distribution is ruled by the behaviour imprinted by 
neutrino decoupling as hot relics at the MeV scale. 
If we denote by P m (y), m being the degree of P m (y), 



P m {y) = Y. c t ) v k 



(4) 



k=0 



the set of polynomials orthonormal with respect to the 
measure y 2 /(exp(y) + 1) 
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it is easy to write the neutrino distribution in terms of 
its moments 
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where we have defined 

m 

^ m = £4 m) Q«T- fc , (7) 

i.e. a linear combination of moments up to order m with 

coefficients ci m . 

For a Fermi-Dirac distribution all moments can be ex- 
pressed in terms of the number density Qa^ or. equiv- 
alently, as functions of the only independent parameter 
T u . The first two moments, namely the neutrino energy 
density defined during the relativistic regime and their 
number density, enter all present analyses of cosmologi- 
cal obscrvables. More precisely, it is customary to define 
the effective number of neutrinos 

\ / a 

where N e g = 3 for a thermal distribution, while it gives 
N c s = 3.04 when including the small effect of neutrino 
heating during the e + — e~ annihilation phase Q- Sim- 
ilarly, Qa are related to the present neutrino energy 
density normalized to the critical density today 

u, u = n„h 2 = 0.058 ^T- 3 QW , (9) 

where m is the sum of the neutrino masses and we have 
assumed for simplicity that the three neutrinos share the 
same distribution. Flavor neutrino oscillations arc in fact 
very effective after BBN, since for temperatures signif- 
icantly lower than McV the background effects of the 
medium are negligible, assuming the neutrino mixing pa- 
rameters favoured by present experimental results, see|4|- 
Therefore, the initial distortions on the flavor neutrino 
spectra are fastly redistributed according to oscillations 
in vacuum. Averaging over all oscillatory terms, it is 
easy to find the final energy spectra of the three neutrino 
mass states. Here we take the approximation of equal 
(non-thermal) distributions for all neutrino mass states. 
In light of present sensitivity of data to neutrino distribu- 
tion, to be discussed in the following, this approximation 
is quite reasonable. 

In the standard (thermal) case we recover the usual 
result 

W " _ 94.1[93.2]eV ' (10) 

where the value between brackets takes into account the 
entropy release to neutrinos from e + — e~ annihilations 
. Actually, there are more parameters which may affect 
the LSS power spectrum, such as the individual neutrino 
masses, but their effect is negligible for a fixed Too- 

What about higher moments? In principle for an ar- 
bitrary distribution function, with the only requirement 
of an exponential decay for large comoving momentum, 
all moments may be relevant in the way this distribu- 
tion is reconstructed by using the expansion ©. In the 



following we will consider only the first two moments 
u) v and N e g as free and independent parameters, to be 
constrained using cosmological data. In fact, we start 
observing that in the limit of a thermal distribution, the 

value of df a (y) is uniquely determined by Qa \ see Eq. 
© . Adopting therefore a conservative approach one may 
think that if no extremely strong deviations from equi- 
librium are present, the expansion © can be truncated 
for some m without losing relevant pieces of information. 
The only way to decide how many parameters should 
be included in the analysis is by studying the sensitiv- 
ity of the available observational data to the distorsion 
of the neutrino distribution. As we will see, using all 
available data from CMB, LSS and possibly BBN, it is 
already very hard to get quite strong constraints on the 
first two moments, so that it appears that presently it 
is not worth comparing more complicated models with 

in) 

observations and introduce further parameters Q a . In 
fact more moments means more parameter degeneracies, 
so that the main conclusion of the analysis would re- 
main unchanged. Of course, in case future observations 
would reach a higher sensitivity on neutrino distribution, 
it would be desirable to include higher order moments, 

(2) 

such as the skewness or the kurtosis, related to Q y a ' and 
Q^a \ respectively, as free parameters in the analysis. 

We now turn to the main point: how the neutrino 
distributions may acquire non-thermal shapes after the 
weak interaction freeze-out. A possible source of such 
non standard contribution may arise from the decays of 
unstable massive particles provided 

i) decays take place out of equilibrium, i.e. for tem- 
peratures smaller than the decaying particle mass 
M . In fact, if decays take place at equilibrium for 
T ~ M, the resulting neutrino distribution would 
be still thermal, possibly with a different ratio of 
neutrino-photon temperatures if the entropy is re- 
leased in different amounts to neutrinos and pho- 
tons. Therefore, in this case the effect consists at 
most in an overall constant factor in neutrino distri- 
bution, changing the energy density in relativistic 
species. 

ii) occur after weak interaction freeze-out. 

The simplest scenario which we will be considering in 
details in the following is the two body decay of a light 
(M < 1 MeV) neutral scalar particle $ whose dynamics 
is dictated by the following lagrangian density 

C= ^d^d^-^MH 2 - ^Y^vm + dnt , (11) 

* i 

with i the flavor index and we have assumed a flavor in- 
dependent coupling A to neutrinos. With £j nt we denote 
the interaction terms of <E> with all other fields except 
neutrinos and responsible for the thermalization of the 
<& quanta at high temperatures. The <E> particles are as- 
sumed to decay in out of equilibrium conditions before 
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4 6 8 10 

Neutrino momentum y 

FIG. 1: Differential number density of relic neutrinos as 
a function of the comoving momentum for the non-thermal 
spectrum in Equation The parameters are A — 0.018, 

y, — 10.5 and a = 1, which corresponds to N e s ~ 4. 



the photon last scattering epoch, so that the produced 
neutrino burst directly influences the CMB anisotropy 
spectrum, as well as the late LSS formation. Different 
scenarios can be envisaged, as for example the case of 
unstable neutrinos fh decaying into a (pseudo) scalar par- 
ticle tp and a lighter neutrino v%, 



(12) 



or three-body decays. It is worth stressing that the aim 
of this paper is to show how combining different cos- 
mologies] observables we can constrain the non-thermal 
contribution to the neutrino background. Despite of the 
fact that we will be mainly dealing for definitcness with 
the out of equilibrium <f> decay case, all our results are 
quite general and can be applied to different scenarios as 
well. Using Eq. I|ll|) we can easily compute the decay 
rate 



T($ -> Vv) 



X 2 M 
8tt 



(13) 



If decays take place out of equilibrium for temperatures 
smaller than the <& mass it follows that the coupling A 
should be very small. In the radiation dominated regime 



A - 10" 



Tp/MeV 
y/AI/MeV 



(14) 



where Tp is the temperature of the thermal part of the 
neutrino background at decay. The coupling A also pro- 
duces scattering processes like <-> $i/ and crossed re- 
actions, which however receive contribution at order A 4 . 
It is easy to see that for the small values implied by Eq. 
(|14f) these processes were never in equilibrium before To, 
so they cannot thermalize the <J> particles with neutrinos 
after decoupling of weak interactions. When the $ par- 
ticles decay, the neutrino distribution gets an additional 



contribution, which in the narrow width limit and in the 
instantaneous decay approximation at Tp corresponds to 
a peaked pulse at = M/(2Td) so that 



y 2 f(y)dy = y 2 



1 



e» + 1 
.4 



dy 



cxp 



y-y* 

2cr 2 



,(15) 



where a -C y*. Relaxing this assumption would be one 
of the many ways to generalize our study. We show an 
example of this non-thermal neutrino spectrum in Figure 
^ where a second peak is clearly seen around y*. With 
this distribution, and accounting for tiny neutrino heat- 
ing effects from e + — e~ annihilations, the lower moments 
expressed in terms of the parameters of Eqs. (JSJ) and (jSJ 
read 



m 



93.2 eV 



2tt 2 

1 + ' 99 3C(3) A 



120 

N eS = 3.04 ( 1 + 0.99— Ay, 



(16) 
(17) 



The value of the constant A, which measures the number 
of neutrinos and antincutrinos produced at decay can be 
expressed in terms of the total number of <f> particles 
before decay in a comoving volume 



A = 



2 n$ 

3 n 



(18) 



with n$ the <E> number density. Since these particles are 
decoupled from the thermal bath since at least the BBN 
epoch, we see that the right hand side of Eq. I|18|) can be 
evaluated at this early epoch, so that the largest value 
of A can be bound by BBN as a function of the $ mass 
and number density. We will discuss this issue in the 
next Section along with the constraints on A and aris- 
ing from present data on BBN, CMB and LSS. Later 
we will show in Section IIVI how future CMB and LSS 
data could improve these constraints. Our results, which 
when expressed in terms of the parameters A and get 
a direct interpretation in the $ particle decay model de- 
scribed so far, nevertheless can be looked upon as quite 
general and give the expected sensitivity that future data 
will likely have in pointing out non-thermal features, if 
any. This sensitivity, as unfortunately it is quite com- 
monly the case, will be crucially depending on the spe- 
cific cosmological model which is assumed as a reference 
one, namely on the number of parameters which enter 
the data fitting procedure. For example, we will stress 
that if we enlarge this set of parameters allowing for ex- 
tra relativistic degrees of freedom in addition to ordinary 
photons and neutrinos, a parameter degeneracy would 
make it very difficult to distinguish between non-thermal 
corrections to the ordinary active neutrino background 
and the presence of extra exotic light particles. 
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III. CURRENT BOUNDS FOR A SCENARIO 
WITH THREE NON-THERMAL NEUTRINOS 

A. Constraints from BBN 

The amount of light nuclei produced during BBN is 
rather sensitive to the value of the Hubble parameter H 
during that epoch, as well as to its time dependence. In 
particular the 4 He mass fraction Y p strongly depends on 
the freeze-out temperature of weak processes which keep 
neutrons and protons in chemical equilibrium. Chang- 
ing the value of H affects the neutron to proton number 
density ratio at the onset of the BBN, which is the key pa- 
rameter entering the final value of Y p and more weakly 
the Deuterium abundance. For a fixed baryon density 
parameter in the range suggested by results of WMAP, 
Lu b = n b ti 2 = 0.023 ±0.001 GJ both nuclei yields are 
in fact monotonically increa sing functions of H . For re- 
cent reviews on BBN see |t! l28i 129^. This fact allows to 
severely constrain any possible extra contribution to the 
total energy density due to other species such as the de- 
coupled massive scalar field $ considered in the present 
analysis. Defining by T$ the corresponding temperature 
parameter, the energy and number density of $ particles 
can be written as 



P<i> 



2^* 



1 



■ dx = 



1/2 



1 



- 1 



C(3) r3 



dx ,(19) 



(20) 



where we have assumed, as from our discussion in the 
previous Section, that the $'s decoupled as hot relics 
well before the onset of BBN. We mention that the effect 
of massive particles with mass of the order of 1 -j- 10 McV 
which are instead coupled to neutrinos or photons during 
BBN has been studied in details in 

H3- The two Eqs. 
(|19f) and 120|) can be combined getting 
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(21) 
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so that the upper bound on p$ from light nuclei abun- 
dances translates into an upper limit to n$ /Tjj and so to 
A, see Eq. as a function of M. To get this bound 

wejase a standard likelihood procedure, as discussed e.g. 
and the BBN numerical code described in details 
|28( . We use the WMAP prior on Wb mentioned before 
and the experimental result of [3lJ for D/H 




D/H = (2.78± 



+0.44n 
0.38; 



xlO" 



(22) 



We notice that, allowing for non-thermal features in neu- 
trino distributions may in principle change the preferred 
value of LOb and the corresponding uncertainty from a 




0.01 



M (MeV) 



FIG. 2: The 1 a (thin lines) and 2 a (thick lines) BBN 
bounds on the $ number density (normalized to T„) versus 
mass M in MeV. The regions above the contours would be 
in disagreement with the observed primordial abundances of 
4 He or D. 



CMB data analysis. However, we will show in the fol- 
lowing that the effect on the baryon density parameter is 
completely negligible and a likelihood analysis provides 
the same value of a standard ACDM model, as consid- 
ered in 01 • It is therefore consistent to adopt this result 
in our BBN study. 

The 4 He mass fraction Y p obtained by extrapolating 
to zero the metallicity measurements performed in dwarf 
irregular and blue compact galaxies is still controversial 
and possibly affected by systematics. There are two dif- 
ferent determinations [33 . l33l | which are only compatible 
by^invoking the large systematic uncertainty quoted in 

Y p = 0.238 ±(0.002) stat ±(0.005) sys , (23) 
Y p = 0.2421 ±(0.0021) stat . (24) 

In the present analysis we consider a conservative esti- 
mate for the experimental 4 He abundance obtained by 
using the results of [34| 



ycxp 
V 



0.245 ±0.007 



(25) 



Finally, we do not use the Li abundance, since the ex- 
perimental estimates of this nuclide may be affected by a 
large depletion mechanism with respect to the primordial 
value, see e.g. [HEl. 

Our results are summarized in Figurc[21 where we show 
the lcr and 2a bounds on n^/T^ obtained using the 4 He 
mass fraction and D/H, assuming that the <&'s decay after 
BBN. As expected the most stringent constraint is due to 
Y p , since this parameter is very sensitive to any change 
of the value of the Hubble parameter. For <!> masses as 
large as 0.1 MeV, the scalars are relativistic during the 
neutron to proton density freeze-out, so we get an al- 
most mass independent bound from Y p of the order of 
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n^jTl < 0.17 (or A < 0.1 at 2a). This corresponds to 
a limit on the non-thermal component of the neutrino 
number density of the order 54 % of the thermal neu- 
trino contribution. On the other hand, for larger masses 
the upper bound on n$ /Tjf decreases approximatively as 
1/M for M larger than few MeV, see Eq. |J2T|). For 
this mass range, the BBN bound from the D abundance 
actually becomes slightly more stringent than that from 
4 He. Note that the BBN bounds still allow a large contri- 
bution of the decay products to the relativistic neutrino 
energy density N g. This value depends quite crucially 
on the ratio M/Td, which measures how far from equi- 
librium conditions the $ decays take place. For values of 
M smaller than 0.1 MeV we get 

#eff< 3.04^1 + 0.1^ , (26) 

where To < 0.01 MeV, since wc assumed that the $ de- 
cays take place after BBN. If To ~ 0.1M or smaller the 
non-thermal component can give a large or even domi- 
nant contribution to N e g. For higher values M > 1 MeV, 
our simple scenario becomes unlikely because it requires 
the <I> decays to take place after BBN, and therefore one 
needs a very large ratio M/Td > 10 2 . However, in this 
case the contribution to iV c ff can still be very large, de- 
spite of the fact that the <& number density rapidly de- 
creases with increasing values of M. In fact, the upper 
bound on N e g scales as 1/Tjj in this regime, sec Eq. (fTTjl . 

B. Adding constraints from cosmological 
perturbations (CMB and LSS) 

We use the public code 2 CAMB [35| in order to com- 
pute the theoretical prediction for the C; coefficients of 
the temperature and polarization power spectra of CMB, 
as well as the matter spectra P{k). We explicitly imple- 
ment into the code the expression of the non-thermal 
phase-space distribution function, in the same fashion 
as massive neutrino chemical potentials were added to 
CMBFAST |3(| in [33]. We focus on three models: the 
minimal ACDM model including three massive neutrinos 
(with equal mass); our extended model with three non- 
thermal massive neutrinos (ACDM+NT) and finally, for 
comparison a model with three ordinary massive ther- 
mal neutrinos plus extra relativistic degrees of freedom 
(ACDM+R), already analyzed in [T3. l3Sj (using a grid- 
based method and a smaller dataset). We compute the 
Bayesian likelihood of each cosmological parameter with 
a Monte Carlo Markov Chain method, using the pub- 
lic code COSMOMC with option MPI on the MAJO- 
rana cluster at Naples University. The likelihood of each 
model is found using three groups of data sets, as im- 
plemented in the public version of COSMOMC: (i) CMB 



2 CAMB Code homepage |http : // camb . info/ 



data: WMAP 0, VSA g(J, CBI 0, ACBAR g2; 
(ii) LSS data: 32 correlated points from 2dFGRS (up 
to fc max = 0.1 h Mpc -1 ) 0131, 14 uncorrected points 
from SDSS g{| (up to the same scale); and (iii) SNIa 
data from Riess et al. g||. For SDSS, we compare di- 
rectly the data with the real space linear matter power 
spectrum multiplied by a scale-independent bias with flat 
prior, which we marginalize out. For 2dF, we include a 
prior on the bias: the linear matter power spectrum (at 
rcdshift z = 0) is rcscalcd by a redshift-space correction 
at z = 0.17 and a bias factor b 

P(k)\7=o.?7 =b 2 (l + \p cS + \p 2 c ^ P{k%ZV P - (27) 

with /? c ff = 0.85/3 and b = n a m 6 / '(3 (see H3 and refer- 
ences therein). For the redshift-space distortion factor 
(5 we adopt the prior (3 = 0.43+0.08 H3. We do not 
include any data from Lyman-a forests 3 . Finally, we im- 
pose the gaussian prior h = 0.72+0.08 (la) from the 
HST Key Project jig, and for the non-thermal model 
we add the BBN prior described in the previous Sec- 
tion: A < 0.1 at the 2-a level (the prior is implemented 
as a half-gaussian centered at A = and with a stan- 
dard deviation of 0.05). As we already said, a remark- 
able feature of the ACDM+NT model is that BBN does 
not impose an upper limit on N e g at the CMB epoch, so 
the results shown for comparison in the ACDM+R case 
will not include an N c g prior, and will not always be 
in agreement with predictions from standard BBN. For 
the ACDM model, the basis of cosmological parameters 
used by the Markov Chain algorithm is as follows: L J 
the overall normalization parameter ln[10 10 7?. r ad] where 
lZ la d is the curvature perturbation in the radiation era; 
2 > the primordial spectrum tilt n s ; 3 ) the baryon den- 
sity lui, = flbh 2 ; 4 ) the (cold+hot) dark matter density 
Wdm = fidm'i 2 ; ^ 0, the ratio of the sound horizon to 
the angular diameter distance multiplied by 100; 6 ' the 
optical depth to reionization r; 7 ) the 2dF redshift-space 
distortion factor (3; 8 - ) the total mass tuq of the three 
neutrinos. In addition, the ACDM+R model includes: 
9 ) the effective neutrino number N c g = 3.04 + AiV re i.. 
Finally, the ACDM+NT model includes: 9) the statisti- 
cal moment Qq 1 ^ of the three non-thermal neutrinos, or 
cquivalently, their relativistic energy density in units of 
that of a standard massless neutrino N e g > 3.04, see Eq. 



3 The non-thermal model includes massive neutrinos which are in- 
dividually more energetic and with a higher average momentum 
than their thermal counterpart. In presence of such a fluid, the 
non-linear small-scale evolution might differ significantly from 
that of ordinary ACDM models. Therefore, we believe that data 
from Lyman-a forests are inappropriate for a conservative anal- 
ysis. We are aware that the bias constraint does rely on obser- 
vations in the mildly non-linear range 0.2 < k < 0.4 h Mpc -1 , 
and we assume that on these scales the non-thermal neutrinos do 
not generate unusually large non-linear corrections which would 
invalidate our bias prior. 
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FIG. 3: Marginalized likelihood for the parameters of each model: ACDM+NT (black, solid), ACDM+R (red, dashed), or 
ACDM (blue, dotted). The independent basis parameters are included in the first ten plots only. The last two plots show the 
related parameters (A, y*) for the non-thermal model. 



©; 10 - ) their statistical moment Qa , or equivalently, 
the dimensionless parameter q = lo u (93.2 eV/mo) > 1, 
see Eqs. © and ifTUj) . The best value of the effective \ 2 
for the three models is shown in the first line of Table [I] 
The ACDM+R and ACDM+NT do not improve signifi- 
cantly the x 2 despite of their larger parameter space (for 
a precise statement on the goodncss-of-fit of the various 
models, one should carry a full "evidence" calculation as 
in 0| ) . On Figure |3 we plot the marginalized likeli- 
hood for each parameter and model. The first ten plots 
correspond to basis (independent) parameters, while the 
last two refer to the derived parameters (A, y*) for the 
non-thermal case, see Eqs. Ill (if and l|17|l . The preferred 
values and l-a errors for each parameter are shown in 
Table D 

For the parameters ln[10 10 7?. ra d], n s , uif,, r, (3, the like- 
lihoods are almost equal in the three cases, fn parti- 
cular, the preferred range for the baryon density in the 
ACDM+NT case is still the one of the standard ACDM 
model, as we already anticipated in Section II. On the 
other hand, in the ACDM+R and ACDM+NT cases, the 
9 distribution is significantly shifted, and Wdm and mo 
have much wider error bars. The existence of a degen- 
eracy direction involving mainly Udm, mo & n d N e g was 



already pointed out in [HEIl for the ACDM+R model, 
and is confirmed here for the two models in which the 
radiation density near the time of decoupling is a free 
parameter 4 . In particular, values as large as Wdm = 0.23, 
mo = 1.5 eV or N g = 9 are still allowed at the 2-er 
level by our set of data and priors (while in the standard 
ACDM case wc find cj dm < 0.130 and m < 0.7 cV at 2- 
cr, in agreement e.g. with Note however that in the 

ACDM+R case, models with iV e ff > 4 are in disagree- 
ment with standard BBN, while in the ACDM+NT case 
they are not. The degeneracy is illustrated by the two- 
dimensional likelihood plots of Figure 0] Finally, we can 
notice that in the ACDM+NT model the bound on the 
parameter q = lj„(93.2 eV/mo) comes essentially from 
the BBN prior A < 0.1, since the CMB and LSS data 
alone would be compatible with much larger deviations 
from a thermal phase-space distribution (up to A = 1 at 
2-er). As a final remark we notice that is poorly con- 



4 This degeneracy is due to the freedom to increase the total radi- 
ation, total matter and critical densities, while keeping fix ed t he 
redshift of equality and tt^. More details can be found in |12| . 
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FIG. 4: Two-dimensional likelihoods illustrating the main de- 
generacy affecting the ACDM+NT model (for the ACDM+R 
model the plots are extremely similar). This degeneracy cor- 
relates N e g mainly with w<im (top left) and mo (middle left), 
and to a smaller extent with 6 (top right) and n s (middle 
right). We also show how N e s is correlated with two related 
parameters: the age of the Universe (bottom left) and the 
Hubble parameter (bottom right). The solid curves refer to 
the 1-ct and 2-a contours of the marginalized likelihood. For 
comparison, the average sample likelihood is shown in color 
levels. 



strained, and large values for it are still allowed. In the 
$ decay scenario this implies that these decays can take 
place in highly out of equilibrium conditions, Tjj <C M, 
the only bound being instead on the scalar particle num- 
ber density at the BBN epoch. 



IV. FUTURE PROSPECTS: CAN WE REMOVE 
THE DEGENERACIES? 

A. Error forecast assuming no extra relativistic 
degrees of freedom 

We will now compute the error expected from future 
experiments on the parameters of the ACDM+NT model. 
Of course, the results that will be derived here are model- 
dependent, since we assume that the cosmological evolu- 
tion can be described with a given set of cosmological pa- 
rameters. The issue of generalization to larger classes of 
models will be considered in the next Section. As usually 
done in this kind of analysis, we first assume that future 





ACDM 


ACDM+R 


ACDM+NT 


Xmin 


1688.2 


1688.0 


1688.0 


ln[10 I,J ft rad ] 


3.2±0.1 


3.2+0.1 


3.2+0.1 


n s 


0.97±0.02 


0.99+0.03 


1.00+0.03 


U>b 


0.0235±0.0010 


0.0231+0.0010 


0.0233+0.0011 




0.121±0.005 


0.17±0.03 


0.17+0.03 


e 


1.043±0.005 


1.033+0.006 


1.033+0.006 


T 


0.13±0.05 


0.13+0.06 


0.15+0.07 


(3 


0.46±0.04 


0.48+0.04 


0.48+0.04 


mo (eV) 


0.3±0.2 


0.8+0.5 


0.7+0.4 


iVefl 


3.04 


6+2 


6+2 


q 


1 


1 


1.25+0.13 


h 


O.o7±0.02 


0.76±0.06 


0.76±0.05 


Age (Gyr) 


13.8±0.2 


12.1+0.9 


12.1+0.8 




0.68±0.03 


0.67+0.03 


0.67+0.03 




14±4 


16+5 


18+6 


£?"8 


0.76±0.06 


0.77+0.07 


0.77+0.07 



TABLE I: Minimum value of the effective \ (defined as 
— 2 In C, where C is the likelihood function) and la confi- 
dence limits for the parameters of the three models under 
consideration. The first ten lines correspond to our basis of 
independent parameters, while the last five refer to related 
parameters. 



data will favor a particular ACDM+NT model that wc 
call the "fiducial model", with cosmological parameter 
values Xi = ir,. Then, we compute the Fisher matrix 



13 dx t dxj 



(28) 



which represents the curvature of the likelihood distri- 
bution C near its maximum for a given set of future ex- 
periments with known sensitivity. Inverting the Fisher 
matrix, we obtain the l-cr error on each parameter, as- 
suming that all other parameters are unknown 



Xi 



(29) 



The present analysis follows the lines of that in [5Cj , and 
we refer the reader to this work for details on the ex- 
pression of the likelihood distribution, for the description 
of future experimental sensitivities and for references of 
seminal papers on the subject. An analysis for massless 
neutrinos and aimed to study the sensitivity on 7V e ff of 
WMAP and Planck experiments has been also discussed 
in |5 1| . When computing numerically the derivative of 
the likelihood with respect to each cosmological param- 
eter, we choose the smallest possible step size above the 
threshold set by the precision of CMBFAST, and we check 
that the final result is almost independent of the exact 
step value. Our first "conservative" set of experiments 
consists in the Planck satellite 5 (without lensing infor- 
mation combined with the completed SDSS redshift 



' |http: //www. rssd. esa. int / index ,php?project=PLAHCK 
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name ln[10 lu ^ rad ] 


n s 








T 


Fhc 


b 2 


m (eV) 




q 




fiducial values 3.1 


0.96 


0.023 


0.14 


0.70 


0.11 


0.24 


1.0 


0.5 


4.0 


l.i 




la error 0.01 


0.009 


0.0003 


0.004 


0.02 


0.005 


0.02 


0.08 


0.3 


0.1 


0.7 


(conservative) 


0.006 


0.003 


0.0001 


0.0007 


0.001 


0.003 


0.004 





0.03 


0.02 


0.08 


(ambitious) 



TABLE II: Expected errors on the parameters of the ACDM+NT model. The first line shows the fiducial values, i.e. the 
parameter values assumed to represent the best fit to the data. The next lines give the associated la errors for the "conservative" 
and "ambitious" sets of future experiments as described in the text. 



survey 6 with effective volume V^ff = 1 h~ 3 Gpc 3 and a 
free bias. We assume that the SDSS data can be com- 
pared with the predictions of linear perturbation the- 
ory up to the scale k niax = 0.1 h Mpc~ x . The second 
"ambitious" set consists in CMBPOL 7 combined with a 
large survey with effective volume V e s = 40 h~ 3 Gpc 3 , no 
bias uncertainty, and fc max = 0.2 h Mpc -1 , which should 
mimic the properties of a large cosmic shear survey like 
the LSST project 8 (this estimate should be regarded 
as indicative only, since for simplicity we compute the 
Fisher matrix from the matter power spectrum instead 
of the spectrum of the projected gravitational potential). 

The results are shown in Table ITT1 for a fiducial model 
comfortably allowed by current bounds, with N e g = 4, 
q = 1.1 and mo = 0.5 eV, corresponding to non-thermal 
correction parameters A = 0.018 and = 10.5. For 
this model, the neutrinos enter the non-relativistic regime 
well after decoupling, and only contribute to 0.5% of the 
critical density today. As a consequence, the two pa- 
rameters q and mo are measurable only from the free- 
streaming effect in the matter power spectrum, while 
all other parameters have a clear signature in the CMB 
anisotropics. More precisely, a combination of these two 
parameters (close to qxmo oc ui v at first approximation) 
controls the amplitude of the free-streaming suppression 
on small scales, while the orthogonal combination con- 
trols the details of the transition on intermediate scales. 
Since LSS experiments do not have enough resolution 
for constraining this last combination, we expect a large 
degeneracy between q and mo, which is confirmed by 
the diagonalization of the Fisher matrix. The direction 
Aq = (Aroo/0.5 eV) appears to be poorly constrained. 
As a consequence, the errors quoted in Table ITT1 for m 
are larger (typically by a factor two) than those expected 
for thermal neutrinos, and q has the biggest relative er- 
ror. Notice that we are not applying in this case any 
prior on q, as was instead the case in the previous Sec- 
tion, where this parameter was related to the $ number 
density at the BBN epoch in the <& to neutrino decay sce- 
nario. Therefore, results in Table |n] are independent of 
any particular model that produces a distortion on the 
relic neutrino spectrum. The remarkable result of this 
Section is that in spite of this degeneracy, future experi- 



http 


//www. sdss 
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//universe 


. nasa . gov/program/inf lat ion . html 
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//www. lsst 
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ments can provide some insight on non-thermal features 
and still measure the neutrino mass with good precision 
even when the thermal assumption is relaxed. The error 
for all parameters but {q, mo) are close to the usual ex- 
pectations based on thermal neutrinos, see for com- 
parison. We also find a sensitivity on 4 He from CMB 
of the same order of magnitude found in |53j| for Planck 
experiment. A small distortion of the first moment de- 
scribed by q ~ 1.16 would be sufficient for a 2-a detec- 
tion with CMBPOL+LSST. Notice that large distortions 
(close to the current upper bounds) would lead to smaller 
errors than those quoted in Table [H] due to their addi- 
tional direct effect on the CMB. 

As a final remark we note that in the CMBPOL+LSST 
scenario the uncertainty on the 4 Hc mass fraction would 
be reduced by a factor two with respect to present sensi- 
tivity from direct measurements of this nuclide via spec- 
troscopic methods. In the framework of the scalar de- 
cay scenario considered in this paper it is worth seeing 
whether this result from CMB may further constrain the 
value of q and/or N c g. The result of course depends on 
the fiducial model which is considered, namely what is 
the value of 4 He which future measurements will suggest 
to be the preferred one. We consider as an example a 
fiducial model predicting Y p = 0.245, in order to compare 
with constraints presently available, see Eq. H25|) . with 
an error of 0.004 (see Table II, "ambitious" set of data). 
Even in this optimistic scenario, comparing this determi- 
nation with the results of the nucleosynthesis theoretical 
code of [2^| we find the bound q < 0.16, a factor two 
larger than what is already obtained by CMB/LSS data 
(0.08, see again Table II). Finally, we recall that in the 
particle decay scenario considered in the paper we cannot 
bound the value of N e g at CMB using BBN, since we are 
assuming that decays take place after BBN and the value 
of N e g can be quite large depending on the M/Td ratio. 
Much stronger constraints come in fact from CMB/LSS 
data alone, order 0.1 and 0.02 in the conservative and 
ambitious scenarios respectively, see Table II. 



B. Degeneracy between non-thermal corrections 
and extra relativistic degrees of freedom 

The bounds derived in the previous Sections depend on 
the assumption that the cosmological observations can 
be described within the framework of the ACDM+NT 
model. In the future, we would be in a position to make 
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such an analysis if we had some independent motiva- 
tions for the introduction of non-thermal distortions in 
the neutrino background. Otherwise, even if the best-fit 
ACDM+NT model turned out to be much more likely 
than the best-fit ACDM, many alternative models would 
compete with each other. 

The ACDM+NT model has a natural "competitor", 
which is the previously discussed ACDM+R model, with 
three standard massive neutrinos and extra relativis- 
tic degrees of freedom. Indeed, in the two models the 
amplitude and scale of the free-streaming effect can be 
tuned independently, which is not usually the case (and 
which cannot be mimicked by the introduction of tilt 
running, spatial curvature, tensor modes, isocurvature 
modes, etc.) 

The results of Section llll Bl suggest that at least within 
current experimental sensitivities, the two models are 
roughly equivalent. However, many theoretical argu- 
ments suggest that the intrinsic differences between the 
two cases can lead to distinct observable signatures. For 
instance, the non-thermal corrections boost the average 
neutrino momentum, so at the level of the equations the 
free-streaming effect is never strictly identical in the two 
cases. 

We have investigated in more details the issue of a pos- 
sible (at least approximate) one-to-one correspondence 
between ACDM+NT and ACDM+R models. We studied 
many examples and found a systematic way to find the 
"twin ACDM+R model" of a given ACDM+NT model: 

1 . for any ACDM+NT model with a small total mass 
too < 0.9 eV, the conclusion is straightforward: 
the ACDM+R model sharing the same values of 
N c g and u> v has identical CMB anisotropy spec- 
tra, and a matter power spectrum differing by less 
than one percent on all scales. Note that the "twin 
ACDM+R model" has a different total neutrino 
mass than the original ACDM+NT one, given by 
the simple relation wig = 9 w[/ VT '' . 

2. for any ACDM+NT model with larger neutrino 
mass - such that the neutrinos are not ultra- 
relativistic at decoupling - the situation is a bit 
more complicated: the "twin ACDM+R model" 
does not share that same values of N e g and ui Vl be- 
cause the time of equality is not strictly the same 
for the two models. Therefore, it is necessary to in- 
crease slightly N c g and J7a in the ACDM+R model 
in order to obtain the same position and amplitude 
of the first peak in the CMB temperature spectrum. 
Finally, a simultaneous increase of N c s, uj m and n s 
makes it possible to match the matter power spec- 
tra without changing the time of equality and the 
shape of the CMB temperature spectrum. 

In both cases the differences are below the percent level 
for the matter power spectrum, even less for the CMB 
anisotropy spectra. This is too small for allowing any 



detection with the next generation of experiments 9 . Of 
course, this conclusion applies to a pure CMB and LSS 
analysis where all parameters are assumed to be un- 
known. Strong evidence for the ACDM+NT could still 
arise from external data sets: for instance, in case of con- 
tradiction (assuming standard neutrinos) between future 
CMB and BBN predictions for the value of N e g; or be- 
tween the neutrino mass scale measured by cosmology 
and by beta decay experiments. 



V. CONCLUSIONS 

In this paper we have studied the possibility that the 
phase space distribution of cosmological relic neutrinos 
may be non-thermal, and how present observations can 
constrain any departure from a standard Fermi-Dirac 
function. The latter is expected to hold to a very good 
accuracy (percent level) in the standard scenario, where 
neutrinos last interact at the freeze out of weak inter- 
actions for temperatures of the order of MeV, and then 
simply free stream keeping an equilibrium distribution. 
However, if neutrino evolution is influenced at some later 
stage by some unknown and exotic phenomena, we expect 
that their distribution could be modified. For example, if 
some decoupled particle eventually decay into neutrinos 
in out of equilibrium conditions and after weak inter- 
action freezing, the neutrino distribution can have large 
departures from a simple Fermi-Dirac function. We have 
investigated this model in details, for the case of light 
(M < 1 MeV) decoupled scalar particles decaying in neu- 
trino/ antineutrino pairs. 

Neutrinos both influence the CMB spectrum and the 
clustering of structures on large scales. Their distribu- 
tion, in fact, enters the value of the matter-radiation 
equality point, as well as for massive neutrinos the mat- 
ter power spectrum suppression at scales smaller than 
the free streaming scale. The neutrino distribution can 
be always fully characterized by the set of its moments 

(n) 

Q a which in principle may all affect the cosmological 
evolution of perturbations. However, present data are 
only sensitive to the first two distribution moments. Ac- 
tually, the total neutrino mass too and the moment Qa ^ 
(corresponding to N c g) are well constrained by the CMB 

and LSS data, while Qa" 1 (corresponding to the ratio 
cu u /mo) is not. This fact, along with the existence of 
a degeneracy among LOdm, N e g and too, implies that it is 
presently impossible to establish whether neutrinos are 
thermally distributed or not, even in the most conser- 
vative scenario of no extra relativistic particles at CMB 
epoch apart from photons and neutrinos. In particular, 



This statement is easy to prove quantitatively by repeating the 
Fisher matrix analysis of the previous section, including an extra 
parameter describing the "linear deformation" of one model into 
another. This parameter results to be far from detectable. 



11 



in the framework of the <& decay scenario described in 
Section [HI the most stringent constraint on Q ® comes 
from BBN only, A < 0.1, rather than from CMB and 
LSS. The non-thermal contribution to the total neutrino 
energy density can then be very large or even dominant, 
since it is ruled by the parameter j/» = Mj (2Tb) which is 
very poorly constrained and can take very large values, 

~ 10, even for A close to the BBN upper bound. 

We have also studied how these results may likely 
change when new data on CMB and LSS will be available. 
Adopting a Fisher matrix analysis we find that combining 
the sensitivity of the Planck experiment with the com- 
pleted SDSS redshift survey, the uncertainty on the neu- 
trino first moment 7V e ff would be reduced by almost an 
order of magnitude, but still the neutrino number den- 
sity will be affected by a large error, see for example the 
results of Table |H] The sensitivity to non-thermal fea- 
tures in the neutrino distribution highly improves only by 
considering a survey with the properties of a shear sur- 
vey like the LSST project, and also the sensitivity of a 
CMB experiment like CMBPOL. Even in this case, how- 
ever, it may be quite involved to test the hypothesis of a 
non thermal neutrino background since, as we stressed in 
Section llVI a purely thermal neutrino background with 
extra relativistic degrees of freedom would provide an 
almost degenerate scenario, with the exception of pre- 
dicting a different energy density at the time of primor- 
dial nucleosynthesis, and a different value for the neu- 
trino mass scale. In the future, evidence for non-thermal 



distorsions could possibly result from an apparent con- 
tradiction between the values of N e g and/or mo probed 
by cosmological perturbations (assuming standard neu- 
trinos) and those indicated respectively by primordial nu- 
cleosynthesis and particle physics experiments. In partic- 
ular, an independent source of information on mg would 
greatly help. Direct neutrino mass measurements like the 
Tritium beta decay experiment KATRIN 1541 or experi- 
ments on neutrinoless double beta decay [53, which are 
expected to reach a sensitivity on the neutrino mass scale 
of the order of 0.1-^0.3 eV, combined with cosmological 
observations could reduce this degeneracy. 

Despite of the fact that by Ockham's razor principle 
it is presently worth assuming a purely thermal distribu- 
tion for relic neutrinos, nevertheless we think that it is re- 
markable that present knowledge of the Ci^B features still 
allows for quite non standard and exotic physics which 
affected the evolution of the Universe at late stages. As 
frequently happened in the past, neutrinos can play a 
role in constraining such possibilities. 
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